%AHA tools for result displaying
%Developed by ruiguo email: ruiguo.ok#163.com

%Read ALL ROIs infomation
%save patiten info to the excel files

% close all ;
clear all ;


T2_cbar = jet;
files_dir = 'H:\Projects\MST1\EXPData\Invivo\AllPatients\MST1';
dirInfo = dir( files_dir );

%% save to excel
% open excel file
str = date ;
T1_T2_w =1;
xlsxfilenm =[ '_Volunteers_AHA_'  str '_T' num2str(T1_T2_w) '.xlsx' ];
filepath = [ files_dir  xlsxfilenm];
try
    Excel=actxgetrunningserver('Excel.application');
catch
    Excel=actxserver('Excel.application');
end

Excel.visible=1;
Workbooks= Excel.Workbooks;
if exist(filepath,'file' )
    Workbook = Excel.Workbooks.Open(filepath);
else
    Workbook = Excel.Workbooks.Add;
    Workbook.SaveAs(filepath);
end
Sheets=Excel.Activeworkbook.Sheets;
Sheet1=Sheets.Item(1);
Sheet1.Activate;

% output all summary info
Start_col = 'A';
row_ix = 1 ;
ActivesheetRange = get(Sheet1,'Range',[(Start_col) num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'SubName');
ActivesheetRange = get(Sheet1,'Range',[Start_col+1 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA1');
ActivesheetRange = get(Sheet1,'Range',[Start_col+2 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA2');
ActivesheetRange = get(Sheet1,'Range',[Start_col+3 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA3');
ActivesheetRange = get(Sheet1,'Range',[Start_col+4 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA4');
ActivesheetRange = get(Sheet1,'Range',[Start_col+5 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA5');
ActivesheetRange = get(Sheet1,'Range',[Start_col+6 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA6');
ActivesheetRange = get(Sheet1,'Range',[Start_col+7 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA7');
ActivesheetRange = get(Sheet1,'Range',[Start_col+8 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA8');
ActivesheetRange = get(Sheet1,'Range',[Start_col+9 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA9');
ActivesheetRange = get(Sheet1,'Range',[Start_col+10 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA10');
ActivesheetRange = get(Sheet1,'Range',[Start_col+11 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA11');
ActivesheetRange = get(Sheet1,'Range',[Start_col+12 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA12');
ActivesheetRange = get(Sheet1,'Range',[Start_col+13 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA13');
ActivesheetRange = get(Sheet1,'Range',[Start_col+14 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA14');
ActivesheetRange = get(Sheet1,'Range',[Start_col+15 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA15');
ActivesheetRange = get(Sheet1,'Range',[Start_col+16 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA16');





for subj_ix = 3:length(dirInfo)
    subject_name = dirInfo(subj_ix).name;
    subjectInfo = dir([files_dir filesep subject_name filesep '*.mat'] );
    image_name =char( regexpi(subjectInfo.name,'.*?(?=.mat)','match'));
    file_dir = [ files_dir filesep subject_name filesep];
    %read ROI info ;
    subject_dir_info = dir( file_dir );
    ROI_state = 0;
    for ix = 3:length(subject_dir_info )
        ROI_dir_name = subject_dir_info(ix).name;
        if(~isempty( regexpi(ROI_dir_name,'(?<=ROIs_).*','match')))
            ROI_dir = [ file_dir ROI_dir_name ];
            ROI_state  = 1;
            break;
        end
    end
    
    if(ROI_state )
        min_slice_num  = 20;
        max_slice_num = 0;
        ROI_dir_info = dir(ROI_dir);
        for ix = 3:length( ROI_dir_info)
            ROI_name = ROI_dir_info(ix).name;
            
            if( ~isempty( regexpi(ROI_name,'S\d\d.mat')) || ~isempty(regexpi(ROI_name,'S\d.mat')))
                str_cell = regexpi(ROI_name,'[\d][\d]{0,}','match') ;
                slice_num = str2num(str_cell{1,1});
                if(min_slice_num>slice_num )
                    min_slice_num=slice_num ;
                end
                if(min_slice_num>slice_num )
                    min_slice_num=slice_num ;
                end
                if(max_slice_num<slice_num )
                    max_slice_num=slice_num ;
                end
            end
        end
        Slice_num_range = [min_slice_num:max_slice_num];
    else
        Slice_num_range = [1:12];
    end
    
    Apical_slice_num = 4;
    Midle_slice_num = 8;
    Base_slice_num = 12;
    Segmentation_tables= {1};
    AHAInfo = {1};
    Thickness_img = [];
    bp_Area = [];
    Start_row =1;  % global setting for the excel output  A D G J
    T1_T2_excel = [ 'T_1';'T_2'];
    T1_T2_str=[ '{T_1}';'{T_2}'];
    T1_T2_DISP = [2000;100];
    % image_name = 'Fit4ParaSlice12sCOSATIRT1T2Fitting_QIAOHUIYU_WIP_3DCOSATIRT1T2_CLEAR_6_1';    % no suffix
    % file_dir =  'D:\paper\CO-register T1 T2 mapping\Results\HealthyVolunteers\PAR REC\Parameters dat\QiaoHuiyu\';

    
    for ll=T1_T2_w:T1_T2_w
        T1_T2 =ll;
      
        Start_row =(subj_ix-2);

        Slice_start_num = Slice_num_range(1,1);
        if(isempty(ROI_dir_name))
            ROI_dir_name = ['ROIs_' image_name];
        end
        for ix = 1 : length(Slice_num_range)
            Obj_slice_num = Slice_num_range(1,ix);
            if(Obj_slice_num>Slice_start_num)
                center_pos(1,1) = mean(Segmentation_tables{1,Slice_start_num}.ENDO_x );
                center_pos(2,1) = mean(Segmentation_tables{1,Slice_start_num}.ENDO_y );
            else
                center_pos = [];
            end
            [rsegment_info,aha_info,Thickness_figure, o_bp_Area]=readsingleroiinfo( file_dir, image_name, Obj_slice_num,center_pos,T1_T2,ROI_dir_name);
            Thickness_figure(Thickness_figure>10000) = 0;
            Thickness_img(:,:,Obj_slice_num) = Thickness_figure;
            bp_Area(:,:,Obj_slice_num) = o_bp_Area;
            Segmentation_tables{1,Obj_slice_num}= rsegment_info;
            AHAInfo{1,Obj_slice_num} = aha_info;
        end
        
        %Generate the AHA data
        % AHA_image = zeros(size(Thickness_figure));
        %% first type
        % AHA_image = Thickness_img(:,:,1) ;
        % Line_st.EPI_x = Segmentation_tables{1,1}.EPI_x+(Segmentation_tables{1,1}.EPI_x-center_pos(1,1))./abs(Segmentation_tables{1,1}.EPI_x-center_pos(1,1));
        % Line_st.EPI_y = Segmentation_tables{1,1}.EPI_y+(Segmentation_tables{1,1}.EPI_y-center_pos(2,1))./abs(Segmentation_tables{1,1}.EPI_y-center_pos(2,1));
        %
        % for ix = 2: length(Slice_num_range)
        %     next_roi_img =  Thickness_img(:,:,ix) ;
        %     next_rsegment_info = Segmentation_tables{1,ix};
        %
        %     pre_roi_img = Thickness_img(:,:,ix-1);
        %     pre_rsegment_info = Segmentation_tables{1,ix-1};
        %
        %     temp_AHA_img = zeros( size( next_roi_img ));
        %     affline_len = length( next_rsegment_info.ENDO_x );
        %     for affline_ix = 1 : affline_len
        %         line_st_x  = Line_st.EPI_x(1,affline_ix);
        %         line_st_y  = Line_st.EPI_y(1,affline_ix);
        %
        %         scale_num = ceil(max( abs(next_rsegment_info.ENDO_x(1,affline_ix)-next_rsegment_info.EPI_x(1,affline_ix)),abs(next_rsegment_info.ENDO_y(1,affline_ix)-next_rsegment_info.EPI_y(1,affline_ix)))/0.5)+1;
        %         if(scale_num > 2 )
        %             affline_x = ceil(linspace(next_rsegment_info.ENDO_x(1,affline_ix),next_rsegment_info.EPI_x(1,affline_ix),scale_num));
        %             affline_y = ceil(linspace(next_rsegment_info.ENDO_y(1,affline_ix),next_rsegment_info.EPI_y(1,affline_ix),scale_num));
        %             for line_ix =  1 : scale_num
        %                 affline_contrast = next_roi_img(affline_y(1,line_ix),affline_x(1,line_ix));
        %                 new_axis_x =  affline_x(1,line_ix) - next_rsegment_info.ENDO_x(1,affline_ix) +line_st_x ;
        %                 new_axis_y =  affline_y(1,line_ix) - next_rsegment_info.ENDO_y(1,affline_ix) +line_st_y;
        %                 temp_AHA_img( ceil(new_axis_y),ceil(new_axis_x)  ) = affline_contrast ;
        %             end
        %         end
        %         Line_st.EPI_x(1,affline_ix) =new_axis_x +abs(new_axis_x)/new_axis_x;
        %         Line_st.EPI_y(1,affline_ix) = new_axis_y +abs(new_axis_y)/new_axis_y;
        %     end
        %     AHA_image(isnan(AHA_image)) = 0 ;
        %     AHA_image = AHA_image +temp_AHA_img;
        % end
        % figure, imagesc( AHA_image );
        
        %% second type
        if(0)
        
        Delta_r = 2;
        St_R = 15 ;
        Total_slice = length(Slice_num_range) ;
        Matrix_R  =  St_R + Total_slice*Delta_r ;
        center_x = Matrix_R +1;
        center_y = Matrix_R +1;
        AHA_matrix = zeros( Matrix_R*2+1,Matrix_R*2+1);
        Sampling_point  = length(Segmentation_tables{1,Slice_start_num}.ENDO_x );
        
        Sample_pt_along_circle = Sampling_point ;
        delta_theta = 2*pi/Sample_pt_along_circle;
        Common_theta_arr =[0:delta_theta:2*pi-delta_theta];
        delta_r = 0.2;
        
        for slice_ix  = 1 :  length(Slice_num_range)
            Obj_slice_num = Slice_num_range(1,slice_ix);
            Common_R_arr = (St_R+(slice_ix-1)*Delta_r):delta_r:(St_R+slice_ix*Delta_r);
            Common_X_axis_arr  = Common_R_arr'*sin(Common_theta_arr)+center_x ;
            Common_Y_axis_arr  = Common_R_arr'*cos(Common_theta_arr)+center_y ;
            affline_len = size(Common_Y_axis_arr,2 );
            next_rsegment_info = Segmentation_tables{1,Obj_slice_num};
            next_roi_img = Thickness_img(:,:,Obj_slice_num);
            temp_AHA_img = zeros( size(AHA_matrix) );
            for affline_ix = 1 : affline_len
                
                scale_num = length( Common_R_arr );
                affline_x = ceil(linspace(next_rsegment_info.ENDO_x(1,affline_ix),next_rsegment_info.EPI_x(1,affline_ix),scale_num));
                affline_y = ceil(linspace(next_rsegment_info.ENDO_y(1,affline_ix),next_rsegment_info.EPI_y(1,affline_ix),scale_num));
                affline_contrast = zeros(1,scale_num);
                for line_ix =  1 : scale_num
                    affline_contrast(1,line_ix) = next_roi_img(affline_y(1,line_ix),affline_x(1,line_ix));
                    %                 new_axis_x =  Common_X_axis_arr(line_ix,affline_ix) ;
                    %                 new_axis_y =  Common_Y_axis_arr(line_ix,affline_ix);
                    %                 temp_AHA_img( fix(new_axis_y),fix(new_axis_x)  ) = affline_contrast ;
                end
                % remove noise
                
                for line_ix =  1 : scale_num
                    if(isnan(affline_contrast(1,line_ix)))
                        affline_contrast(1,line_ix) = mean( affline_contrast(~isnan(affline_contrast)) );
                    end
                    %                 new_axis_x =  Common_X_axis_arr(line_ix,affline_ix) ;
                    %                 new_axis_y =  Common_Y_axis_arr(line_ix,affline_ix);
                    %                 temp_AHA_img( fix(new_axis_y),fix(new_axis_x)  ) = affline_contrast ;
                end
                for line_ix =  1 : scale_num
                    %                 affline_contrast(1,line_ix) = next_roi_img(affline_y(1,line_ix),affline_x(1,line_ix));
                    new_axis_x =  Common_X_axis_arr(line_ix,affline_ix) ;
                    new_axis_y =  Common_Y_axis_arr(line_ix,affline_ix);
                    temp_AHA_img( fix(new_axis_y),fix(new_axis_x)  ) =  mean( affline_contrast(~isnan(affline_contrast)) ) ;
                end
            end
            BIN_AHA_matrix = zeros( size(AHA_matrix));
            BIN_AHA_matrix(AHA_matrix == 0) = 1;
            BIN_temp_AHA_matrix = zeros( size(AHA_matrix));
            BIN_temp_AHA_matrix(temp_AHA_img > 0) = 1;
            Mask = BIN_AHA_matrix.*BIN_temp_AHA_matrix ;
            AHA_matrix = AHA_matrix + temp_AHA_img.*Mask ;
        end
        
        
        f = figure('Visible','Off');,
        hold on;
        axis equal;
        set(gca, 'Color','white');
        imagesc(AHA_matrix,[10 T1_T2_DISP(T1_T2)]), colormap( T2_cbar );
        %plot 45D line
        R_apical_st = St_R ;
        R_apical_ed = Delta_r*Apical_slice_num+R_apical_st;
        Common_R_arr =R_apical_st:Delta_r: R_apical_ed;
        Common_theta_arr  = pi/4:pi/2:2*pi;
        Common_X_axis_arr  = Common_R_arr'*sin(Common_theta_arr)+center_x ;
        Common_Y_axis_arr  = Common_R_arr'*cos(Common_theta_arr)+center_y ;
        % for pt_line_ix = 1:size(Common_X_axis_arr,2)
        %     plot(Common_X_axis_arr(:,pt_line_ix),Common_Y_axis_arr(:,pt_line_ix),'k','LineWidth',3 );
        % end
        % plot circle
        
        R_mild_st =  Delta_r*Apical_slice_num+R_apical_st ;
        R_mild_ed = Delta_r*Total_slice+St_R;
        Common_R_arr =R_mild_st:Delta_r: R_mild_ed;
        Common_theta_arr  = pi/2:pi/2:2.6*pi;
        Common_X_axis_arr  = Common_R_arr'*sin(Common_theta_arr)+center_x ;
        Common_Y_axis_arr  = Common_R_arr'*cos(Common_theta_arr)+center_y ;
        % for pt_line_ix = 1:size(Common_X_axis_arr,2)
        %     plot(Common_X_axis_arr(:,pt_line_ix),Common_Y_axis_arr(:,pt_line_ix),'k','LineWidth',3 );
        % end
        %plot start
        
        % Common_theta_arr  = linspace(0,2*pi, 100);
        % R_midle_edge = St_R ;
        % Common_X_axis_arr  = R_midle_edge*sin(Common_theta_arr)+center_x ;
        % Common_Y_axis_arr  = R_midle_edge*cos(Common_theta_arr)+center_y ;
        %  plot(Common_X_axis_arr,Common_Y_axis_arr,'k','LineWidth',3 );
        %
        % %plot apical
        % Common_theta_arr  = linspace(0,2*pi, 100);
        % R_apical_edge = R_apical_ed;
        % Common_X_axis_arr  = R_apical_edge*sin(Common_theta_arr)+center_x ;
        % Common_Y_axis_arr  = R_apical_edge*cos(Common_theta_arr)+center_y ;
        %  plot(Common_X_axis_arr,Common_Y_axis_arr,'k','LineWidth',3 );
        %
        % %plot midle
        % Common_theta_arr  = linspace(0,2*pi, 100);
        % R_midle_edge = Midle_slice_num*Delta_r +St_R ;
        % Common_X_axis_arr  = R_midle_edge*sin(Common_theta_arr)+center_x ;
        % Common_Y_axis_arr  = R_midle_edge*cos(Common_theta_arr)+center_y ;
        %  plot(Common_X_axis_arr,Common_Y_axis_arr,'k','LineWidth',3 );
        %
        % %plot base
        %
        % Common_theta_arr  = linspace(0,2*pi, 100);
        % R_midle_edge = Base_slice_num*Delta_r +St_R ;
        % Common_X_axis_arr  = R_midle_edge*sin(Common_theta_arr)+center_x ;
        % Common_Y_axis_arr  = R_midle_edge*cos(Common_theta_arr)+center_y ;
        %  plot(Common_X_axis_arr,Common_Y_axis_arr,'k','LineWidth',3 );
        axis off;
        hold off;
        print(f,'-dpng','-r300', [ file_dir filesep  ['AHA' T1_T2_str(T1_T2,:)] '.png']);
        print(f,'-deps','-r300', [ file_dir filesep  ['AHA' T1_T2_str(T1_T2,:)] '.eps']);
        %% third type, slice result
        roi_image = zeros( size(Thickness_img,1 ),size(Thickness_img,2));
        Mean_std_vs = zeros( 2, Slice_num_range(end) );
        Sub_Large_zero = [];
        bp_whole_LV = [];
        for ix = 1 : length(Slice_num_range)
            Obj_slice_num = Slice_num_range(1,ix);
            roi_image =Thickness_img(:,:,Obj_slice_num);
            bp_Img =  bp_Area(:,:,Obj_slice_num);
            bp_Img_noZero = bp_Img(bp_Img>0 &  bp_Img<5000);
            Large_zero = roi_image(roi_image>0 &  roi_image<3000 );
            Mean_v = mean( Large_zero );
            Std_v = std( Large_zero );
            %     Large_zero(Large_zero > (  Mean_v+3*Std_v )) = 0 ;
            %     Large_zero(Large_zero < (  Mean_v-3*Std_v )) = 0 ;
            sub_Large_zero = Large_zero(Large_zero>0);
            Sub_Large_zero(1,end+1:end+length(sub_Large_zero)) = sub_Large_zero ;
            bp_whole_LV(1,end+1:end+length(bp_Img_noZero)) = bp_Img_noZero; 
            
            Mean_std_vs(1,Obj_slice_num ) = mean( sub_Large_zero );
            Mean_std_vs(2,Obj_slice_num ) = std( sub_Large_zero );
        end
        f = figure('Visible','Off');,
        % title('Myocardium T1 over slices( Mean �� Std )','fontsize',20);
        hold on;
        box on;
        for ix = 1 : length(Slice_num_range)
            Obj_slice_num = Slice_num_range(1,ix);
            mean_v  = Mean_std_vs(1,Obj_slice_num);
            std_v  = Mean_std_vs(2,Obj_slice_num);
            plot([Obj_slice_num-0.2 Obj_slice_num+0.2],[mean_v-std_v mean_v-std_v ],'k','linewidth',3);
            plot([Obj_slice_num Obj_slice_num],[mean_v-std_v mean_v+std_v ],'k','linewidth',3);
            plot([Obj_slice_num-0.2 Obj_slice_num+0.2],[mean_v+std_v mean_v+std_v ],'k','linewidth',3);
        end
        plot(Slice_num_range,Mean_std_vs(1,Slice_start_num:end),'k','linewidth',3);
        xlabel('Slice number','FontName','Calibrit','fontsize',20);
        ylabel( [T1_T2_str(T1_T2,:) '[ms]'],'FontName','Calibri','fontsize',20);
        set(gca,'FontName','Calibri', 'FontSize', 16);
        set(gca,'XTick',[0 Slice_num_range]);
        set(gca, 'LineWidth',2);
        axis([Slice_num_range(1,1)-1,Slice_num_range(1,end)+1,min(Mean_std_vs(1,Slice_start_num:end))-max(Mean_std_vs(2,Slice_start_num:end)),max(Mean_std_vs(1,Slice_start_num:end))+max(Mean_std_vs(2,Slice_start_num:end))]);
        hold off;
        print(f,'-dpng','-r300', [ file_dir filesep  ['T1_T2 over slices' T1_T2_str(T1_T2,:)] '.png']);
        print(f,'-deps','-r300', [ file_dir filesep  ['T1_T2 over slices' T1_T2_str(T1_T2,:)] '.eps']);
        %guass plot
        if(T1_T2==1)
            his_gap =10;
        else
            his_gap = 1;
        end
        T1roi = Sub_Large_zero;
        x = min(T1roi(find(T1roi>0))):his_gap:max(T1roi(find(T1roi>0)));
        x = x';
        h = histc(T1roi(find(T1roi>0)),x); % centered on integers
        hs = h;
        [sigma,mu,A] = customGaussFit(x,hs);
        f = figure('Visible','Off');,
        hold on;
        box on;
        bar(x,hs,1) % centered on integers
        mu = mean(T1roi);
        left_leaf =abs( min(T1roi(find(T1roi>0))) -  mu);
        right_leaf =abs( max(T1roi(find(T1roi>0))) -  mu);
        max_x_axis = max( left_leaf,right_leaf );
        sigma = std(T1roi);
        x = ( mu -max_x_axis ):( mu + max_x_axis );
        plot(x, A*exp(-(x-mu).^2/(2*sigma^2)),'--r','linewidth',2)
        titlename = ['Total: ', num2str(sum(h)), ' pixels, Average ', num2str(round(mu)),...
            ' ms, \sigma = ', num2str(round(sigma))]
        title(titlename)
        ylabel('Number of pixels','FontName','Calibri','fontsize',20)
        xlabel( [ T1_T2_str(T1_T2,:) '[ms]' ],'FontName','Calibri','fontsize',20);
        set(gca,'FontName','Calibri', 'FontSize', 16);
        set(gca, 'LineWidth',2);
        hold off;
        print(f,'-dpng','-r300', [ file_dir filesep  ['guass' T1_T2_str(T1_T2,:)] 'mean' num2str(round(mu)) 'sigma' num2str(round(sigma))  '.png']);
        print(f,'-deps','-r300', [ file_dir filesep  ['guass' T1_T2_str(T1_T2,:)] 'mean' num2str(round(mu)) 'sigma' num2str(round(sigma))  '.eps']);
        
        end
        
        %% for AHA statisitics
        AHA_seg =1;
        if(AHA_seg == 1 )
            AHAseg_statistics = {};
            Apex_end_slice_num = 1;
            Four_seg_cnt = 0;
            Six_seg_cnt = 0;
            
            for ix = 1 : length(Slice_num_range)
                Obj_slice_num = Slice_num_range(1,ix);
                %get ROI  image
                ROI_figure = Thickness_img(:,:,Obj_slice_num) ;
                %get ROI info
                segment_info = Segmentation_tables{1,Obj_slice_num};
                %AHA indo
                aha_info = AHAInfo{1,Obj_slice_num};
                
                center_x = mean( segment_info.ENDO_x );
                center_y = mean( segment_info.ENDO_y );
                
                EPI_R = sqrt( ( segment_info.EPI_x -center_x).^2+( segment_info.EPI_y -center_y ).^2);
                Max_EPI_R = max(EPI_R);
                Max_EPI_R= Max_EPI_R(1,1);
                SEG_angle = -2*pi/aha_info.seg;
                Common_x = [0:Max_EPI_R];
                Common_y = zeros(size( Common_x ));
                aha_seg_avg=zeros(aha_info.seg,1);
                aha_seg_std=zeros(aha_info.seg,1);
                
                
                seg_area = {};
                for seg_i = 1:aha_info.seg
                    
                    angle1 = (seg_i - 2)*SEG_angle+ aha_info.angle;
                    angle2 = (seg_i - 1)*SEG_angle+ aha_info.angle;
        
                    [Line_1_x, Line_1_y]= rotatepoint( Common_x ,Common_y, angle1 );
                    [Line_2_x, Line_2_y]= rotatepoint( Common_x ,Common_y, angle2 );
                    arc_angle = linspace(angle1,angle2,10);
                    arc_angle_x = zeros(length( arc_angle)+1,1);
                    arc_angle_y = zeros(length( arc_angle)+1,1);
                    for arc_i = 1:length(arc_angle)
                        [temp_line_x,temp_line_y ] = rotatepoint( Common_x ,Common_y, arc_angle(arc_i) );
                        arc_angle_x( arc_i ) = temp_line_x(end)+center_x;
                        arc_angle_y( arc_i ) = temp_line_y(end)+center_y;
                    end
                    arc_angle_x(end) =  center_x ;
                    arc_angle_y(end) =  center_y ;
                    %ploy the arc
                    blankImg = zeros(size(ROI_figure,1),size(ROI_figure,2));
                    [mesh_x,mesh_y] = meshgrid([1:size(blankImg,2)],[1:size(blankImg,1)]);
                    arc_RR = inpolygon(mesh_x,mesh_y,arc_angle_x,arc_angle_y);
                    [arc_i,arc_j] = find(arc_RR);
                    blankImg(sub2ind(size(blankImg),arc_i,arc_j)) = 1;
                    blankImg = blankImg.*ROI_figure;
                    %         aha_seg_avg(seg_i)=mean(blankImg( blankImg > 0 & blankImg < 5000 ));
                    %         aha_seg_std(seg_i)=std(blankImg( blankImg > 0 & blankImg < 5000 )) ;
                    
                    aha_seg_avg(seg_i)=mean(blankImg( blankImg > 0  ));
                    aha_seg_std(seg_i)=std(blankImg( blankImg > 0     )) ;
                    seg_area{1,seg_i} = blankImg( blankImg > 0  );
                    %         figure,
                    %         imagesc(T1(:,:,Obj_slice_num,2),[ 0 100]);
                    %         hold on;
                    %           plot(Line_1_x+center_x,Line_1_y+center_y,'r-.')
                    %            plot(Line_2_x+center_x,Line_2_y+center_y,'g-.')
                    %         imagesc( ROI_figure );
                    %         plot(arc_j,arc_i,'g+' );
                    %         imagesc(blankImg);
                    %         hold off;
                end
                % recode the apex pos
                if( aha_info.seg == 4 )
                    Apex_end_slice_num =  Obj_slice_num ;
                    Four_seg_cnt = Four_seg_cnt + 1;
                    slice_info.settal_avg = mean([ seg_area{1,2}]);
                    slice_info.settal_std = std([ seg_area{1,2}]);
                elseif( aha_info.seg == 6 )
                    slice_info.settal_avg = mean([ seg_area{1,2}'  seg_area{1,3}']);
                    slice_info.settal_std = std([ seg_area{1,2}'  seg_area{1,3}']);
                    Six_seg_cnt = Six_seg_cnt + 1;
                end
                slice_info.aha_info = aha_info;
                slice_info.aha_seg_avg = aha_seg_avg;
                slice_info.aha_seg_std = aha_seg_std;
                slice_info.area = seg_area;
                AHAseg_statistics{1,Obj_slice_num} = slice_info ;
            end
            % for base statisitcs
            AHA16seg = {1,16};
            Area_13 = [];
            Area_14 = [];
            Area_15 = [];
            Area_16 = [];
            for ix = 1 : Four_seg_cnt
                Obj_slice_num = Slice_num_range(1,ix);
                slice_info = AHAseg_statistics{1,Obj_slice_num};
                Area_13 = [ Area_13 slice_info.area{1,1}' ];
                Area_14 = [ Area_14 slice_info.area{1,2}' ];
                Area_15 = [ Area_15 slice_info.area{1,3}' ];
                Area_16 = [ Area_16 slice_info.area{1,4}' ];
            end
            %3sigma principle
            % Area_13(Area_13>( mean(Area_13)+3*std(Area_13)) | Area_13<( mean(Area_13)-3*std(Area_13))) = 0;
            % Area_14(Area_14>( mean(Area_14)+3*std(Area_14)) | Area_14<( mean(Area_14)-3*std(Area_14))) = 0;
            % Area_15(Area_15>( mean(Area_15)+3*std(Area_15)) | Area_15<( mean(Area_15)-3*std(Area_15))) = 0;
            % Area_16(Area_16>( mean(Area_16)+3*std(Area_16)) | Area_16<( mean(Area_16)-3*std(Area_16))) = 0;
            
            seg.mean = mean(Area_13());
            seg.std= std(Area_13);
            AHA16seg{1,13} = seg;
            seg.mean = mean(Area_14);
            seg.std= std(Area_14);
            AHA16seg{1,14} = seg;
            seg.mean = mean(Area_15);
            seg.std= std(Area_15);
            AHA16seg{1,15} = seg;
            seg.mean = mean(Area_16);
            seg.std= std(Area_16);
            AHA16seg{1,16} = seg;
            
            % for midle statisitcs
            Area_7 = [];
            Area_8 = [];
            Area_9 = [];
            Area_10 = [];
            Area_11 = [];
            Area_12 = [];
            Middle_slice = fix(( length(Slice_num_range)-Four_seg_cnt )/2);
            for ix = Four_seg_cnt+1:Middle_slice+Four_seg_cnt
                Obj_slice_num = Slice_num_range(1,ix);
                slice_info = AHAseg_statistics{1,Obj_slice_num};
                Area_7 = [ Area_7 slice_info.area{1,1}' ];
                Area_8 = [ Area_8 slice_info.area{1,2}' ];
                Area_9= [ Area_9 slice_info.area{1,3}' ];
                Area_10 = [ Area_10 slice_info.area{1,4}' ];
                Area_11 = [ Area_11 slice_info.area{1,5}' ];
                Area_12 = [ Area_12 slice_info.area{1,6}' ];
            end
            % 3 times sigma principle
            % Area_12(Area_12>( mean(Area_12)+3*std(Area_12)) | Area_12<( mean(Area_12)-3*std(Area_12))) = 0;
            % Area_11(Area_11>( mean(Area_11)+3*std(Area_11)) | Area_11<( mean(Area_11)-3*std(Area_11))) = 0;
            % Area_10(Area_10>( mean(Area_10)+3*std(Area_10)) | Area_10<( mean(Area_10)-3*std(Area_10))) = 0;
            % Area_9(Area_9>( mean(Area_9)+3*std(Area_9)) | Area_9<( mean(Area_9)-3*std(Area_9))) = 0;
            % Area_8(Area_8>( mean(Area_8)+3*std(Area_8)) | Area_8<( mean(Area_8)-3*std(Area_8))) = 0;
            % Area_7(Area_7>( mean(Area_7)+3*std(Area_7)) | Area_7<( mean(Area_7)-3*std(Area_7))) = 0;
            
            seg.mean = mean(Area_12);
            seg.std= std(Area_12);
            AHA16seg{1,12} = seg;
            seg.mean = mean(Area_11);
            seg.std= std(Area_11);
            AHA16seg{1,11} = seg;
            seg.mean = mean(Area_10);
            seg.std= std(Area_10);
            AHA16seg{1,10} = seg;
            seg.mean = mean(Area_9);
            seg.std= std(Area_9);
            AHA16seg{1,9} = seg;
            seg.mean = mean(Area_8);
            seg.std= std(Area_8);
            AHA16seg{1,8} = seg;
            seg.mean = mean(Area_7);
            seg.std= std(Area_7);
            AHA16seg{1,7} = seg;
            % for midle statisitcs
            Area_1 = [];
            Area_2 = [];
            Area_3 = [];
            Area_4 = [];
            Area_5 = [];
            Area_6 = [];
            for ix = Middle_slice+Four_seg_cnt+1:length(Slice_num_range)
                Obj_slice_num = Slice_num_range(1,ix);
                slice_info = AHAseg_statistics{1,Obj_slice_num};
                Area_1 = [ Area_1 slice_info.area{1,1}' ];
                Area_2 = [ Area_2 slice_info.area{1,2}' ];
                Area_3= [ Area_3 slice_info.area{1,3}' ];
                Area_4 = [ Area_4 slice_info.area{1,4}' ];
                Area_5 = [ Area_5 slice_info.area{1,5}' ];
                Area_6 = [ Area_6 slice_info.area{1,6}' ];
            end
            % Area_6(Area_6>( mean(Area_6)+3*std(Area_6)) | Area_6<( mean(Area_6)-3*std(Area_6))) = 0;
            % Area_5(Area_5>( mean(Area_5)+3*std(Area_5)) | Area_5<( mean(Area_5)-3*std(Area_5))) = 0;
            % Area_4(Area_4>( mean(Area_4)+3*std(Area_4)) | Area_4<( mean(Area_4)-3*std(Area_4))) = 0;
            % Area_3(Area_3>( mean(Area_3)+3*std(Area_3)) | Area_3<( mean(Area_3)-3*std(Area_3))) = 0;
            % Area_2(Area_2>( mean(Area_2)+3*std(Area_2)) | Area_2<( mean(Area_2)-3*std(Area_2))) = 0;
            % Area_1(Area_1>( mean(Area_1)+3*std(Area_1)) | Area_1<( mean(Area_1)-3*std(Area_1))) = 0;
            
            seg.mean = mean(Area_6);
            seg.std= std(Area_6);
            AHA16seg{1,6} = seg;
            seg.mean = mean(Area_5);
            seg.std= std(Area_5);
            AHA16seg{1,5} = seg;
            seg.mean = mean(Area_4);
            seg.std= std(Area_4);
            AHA16seg{1,4} = seg;
            seg.mean = mean(Area_3);
            seg.std= std(Area_3);
            AHA16seg{1,3} = seg;
            seg.mean = mean(Area_2);
            seg.std= std(Area_2);
            AHA16seg{1,2} = seg;
            seg.mean = mean(Area_1);
            seg.std= std(Area_1);
            AHA16seg{1,1} = seg;
        end
        
        
        %export to the excel
        % start_row = 1 ;
%         row_ix =  Start_row ;
%         Start_col = 'A';
%         Sub_item = 'B';
%         Dat_st_col = 'C';
%         Merage_row = 1;
%         
%         Merage_info = [ Start_col num2str(row_ix) ':' char(Start_col+10) num2str(row_ix)];
%         Sheet1.Range(Merage_info).MergeCells=1;
%         ActivesheetRange = get(Sheet1,'Range',Merage_info);
%         set(ActivesheetRange, 'Value', subject_name);
%         row_ix = row_ix+1 ;
%         %save name
%         Merage_info = [ Start_col num2str(row_ix) ':' char(Start_col+10) num2str(row_ix)];
%         Sheet1.Range(Merage_info).MergeCells=1;
%         ActivesheetRange = get(Sheet1,'Range',Merage_info);
%         set(ActivesheetRange, 'Value', T1_T2_excel(T1_T2,:));
%         
%         row_ix = row_ix+1 ;
%         ActivesheetRange = get(Sheet1,'Range',[(Start_col) num2str(row_ix) ]);
%         set(ActivesheetRange, 'Value', 'Slice num');
%         ActivesheetRange = get(Sheet1,'Range',[Start_col num2str(row_ix+1) ]);
%         set(ActivesheetRange, 'Value', 'Avg');
%         ActivesheetRange = get(Sheet1,'Range',[Start_col num2str(row_ix+2) ]);
%         set(ActivesheetRange, 'Value', 'SD');
%         ActivesheetRange = get(Sheet1,'Range',[Start_col num2str(row_ix+3) ]);
%         set(ActivesheetRange, 'Value', 'Septal_avg');
%         ActivesheetRange = get(Sheet1,'Range',[Start_col num2str(row_ix+4) ]);
%         set(ActivesheetRange, 'Value', 'Septal_avg_std');
%         
%         ActivesheetRange = get(Sheet1,'Range',[Start_col num2str(row_ix+5) ]);
%         set(ActivesheetRange, 'Value', 'Septal_avg');
%         ActivesheetRange = get(Sheet1,'Range',[Start_col num2str(row_ix+6) ]);
%         set(ActivesheetRange, 'Value', 'Septal_avg_std');
%         
%         
%         septal_avg = [];
%         septal_std = [];
%         for ix = 1 : length(Slice_num_range)
%             Obj_slice_num = Slice_num_range(1,ix);
%             ActivesheetRange = get(Sheet1,'Range',[char(Start_col+ix) num2str(row_ix) ]);
%             set(ActivesheetRange, 'Value',Obj_slice_num);
%             mean_v  = Mean_std_vs(1,Obj_slice_num);
%             ActivesheetRange = get(Sheet1,'Range',[char(Start_col+ix) num2str(row_ix+1) ]);
%             set(ActivesheetRange, 'Value', mean_v);
%             std_v  = Mean_std_vs(2,Obj_slice_num);
%             ActivesheetRange = get(Sheet1,'Range',[char(Start_col+ix) num2str(row_ix+2) ]);
%             set(ActivesheetRange, 'Value', std_v );
%             if(AHA_seg)
%                 septal_avg = [ septal_avg AHAseg_statistics{1,Obj_slice_num}.settal_avg];
%                 ActivesheetRange = get(Sheet1,'Range',[char(Start_col+ix) num2str(row_ix+3) ]);
%                 set(ActivesheetRange, 'Value', septal_avg(ix) );
%                 septal_std =[septal_std  AHAseg_statistics{1,Obj_slice_num}.settal_std];
%                 ActivesheetRange = get(Sheet1,'Range',[char(Start_col+ix) num2str(row_ix+4) ]);
%                 set(ActivesheetRange, 'Value', septal_std(ix) );
%             end
%         end
        
%         ActivesheetRange = get(Sheet1,'Range',[char(Start_col+length(Slice_num_range)+1 ) num2str(row_ix) ]);
%         set(ActivesheetRange, 'Value','AVG');
%         ActivesheetRange = get(Sheet1,'Range',[char(Start_col+length(Slice_num_range)+1) num2str(row_ix+1) ]);
%         set(ActivesheetRange, 'Value', mean( Mean_std_vs(1,Slice_num_range)));
%         
%         ActivesheetRange = get(Sheet1,'Range',[char(Start_col+length(Slice_num_range)+1) num2str(row_ix+2) ]);
%         set(ActivesheetRange, 'Value', std( Mean_std_vs(1,Slice_num_range)));
        
        if(0)
            ActivesheetRange = get(Sheet1,'Range',[char(Start_col+length(Slice_num_range)+1) num2str(row_ix+3) ]);
            set(ActivesheetRange, 'Value', mean(  septal_avg));
            
            ActivesheetRange = get(Sheet1,'Range',[char(Start_col+length(Slice_num_range)+1) num2str(row_ix+4) ]);
            set(ActivesheetRange, 'Value', std( septal_avg));    %% for std
            
            ActivesheetRange = get(Sheet1,'Range',[char(Start_col+length(Slice_num_range)+2) num2str(row_ix+3) ]);
            set(ActivesheetRange, 'Value', mean(  [Area_8, Area_9]  ));
            
            ActivesheetRange = get(Sheet1,'Range',[char(Start_col+length(Slice_num_range)+2) num2str(row_ix+4) ]);
            set(ActivesheetRange, 'Value', std( [Area_8, Area_9] ));    %% for std
            
        end
       
        
        
        if(0)
            row_ix = row_ix+7;
            Merage_info = [ Start_col num2str(row_ix) ':' char(Start_col+10) num2str(row_ix)];
            Sheet1.Range(Merage_info).MergeCells=1;
            ActivesheetRange = get(Sheet1,'Range',Merage_info);
            set(ActivesheetRange, 'Value', [T1_T2_excel(T1_T2,:)  '16AHAsegments']);
            ActivesheetRange = get(Sheet1,'Range',[(Start_col) num2str(row_ix+1) ]);
            set(ActivesheetRange, 'Value', 'AHA num');
            ActivesheetRange = get(Sheet1,'Range',[Start_col num2str(row_ix+2) ]);
            set(ActivesheetRange, 'Value', 'Avg');
            ActivesheetRange = get(Sheet1,'Range',[Start_col num2str(row_ix+3) ]);
            set(ActivesheetRange, 'Value', 'SD');
            
            for ix = 1 : 16
                seg = AHA16seg{1,ix};
                ActivesheetRange = get(Sheet1,'Range',[char(Start_col+ix) num2str(row_ix+1) ]);
                set(ActivesheetRange, 'Value',ix);
                
                ActivesheetRange = get(Sheet1,'Range',[char(Start_col+ix) num2str(row_ix+2) ]);
                set(ActivesheetRange, 'Value', seg.mean);
                
                ActivesheetRange = get(Sheet1,'Range',[char(Start_col+ix) num2str(row_ix+3) ]);
                set(ActivesheetRange, 'Value', seg.std);
            end
        end
    end
    

%export summary info   
Start_col = 'A';
row_ix =subj_ix-2;

ActivesheetRange = get(Sheet1,'Range',[(Start_col) num2str(row_ix+1) ]);
set(ActivesheetRange, 'Value', subject_name);

 if(AHA_seg)
    for ix = 1 : 16 
    seg = AHA16seg{1,ix};    
    ActivesheetRange = get(Sheet1,'Range',[char(Start_col+ix) num2str(row_ix+1) ]);
    set(ActivesheetRange, 'Value', seg.mean);     
    end
 end

Start_col = 'A';
row_ix =subj_ix-2+length(dirInfo)+3;

ActivesheetRange = get(Sheet1,'Range',[(Start_col) num2str(row_ix+1) ]);
set(ActivesheetRange, 'Value', subject_name);

 if(AHA_seg)
    for ix = 1 : 16 
    seg = AHA16seg{1,ix};    
    ActivesheetRange = get(Sheet1,'Range',[char(Start_col+ix) num2str(row_ix+1) ]);
    set(ActivesheetRange, 'Value', seg.std);     
    end
 end
 
end

Start_col = 'A';
row_ix = length(dirInfo)+4 ;
ActivesheetRange = get(Sheet1,'Range',[(Start_col) num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'SubName');
ActivesheetRange = get(Sheet1,'Range',[Start_col+1 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA1');
ActivesheetRange = get(Sheet1,'Range',[Start_col+2 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA2');
ActivesheetRange = get(Sheet1,'Range',[Start_col+3 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA3');
ActivesheetRange = get(Sheet1,'Range',[Start_col+4 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA4');
ActivesheetRange = get(Sheet1,'Range',[Start_col+5 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA5');
ActivesheetRange = get(Sheet1,'Range',[Start_col+6 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA6');
ActivesheetRange = get(Sheet1,'Range',[Start_col+7 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA7');
ActivesheetRange = get(Sheet1,'Range',[Start_col+8 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA8');
ActivesheetRange = get(Sheet1,'Range',[Start_col+9 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA9');
ActivesheetRange = get(Sheet1,'Range',[Start_col+10 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA10');
ActivesheetRange = get(Sheet1,'Range',[Start_col+11 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA11');
ActivesheetRange = get(Sheet1,'Range',[Start_col+12 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA12');
ActivesheetRange = get(Sheet1,'Range',[Start_col+13 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA13');
ActivesheetRange = get(Sheet1,'Range',[Start_col+14 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA14');
ActivesheetRange = get(Sheet1,'Range',[Start_col+15 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA15');
ActivesheetRange = get(Sheet1,'Range',[Start_col+16 num2str(row_ix) ]);
set(ActivesheetRange, 'Value', 'AHA16');

%save statical results
Workbook.SaveAs(filepath);
invoke(Excel, 'Quit');
delete(Excel);


